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Abstract 

We present a fast algorithm for generating full sky, high resolution (~ 5') simulations of 
the CMB anisotropy pattern. We also discuss the inverse problem, that of evaluating from 
such a map the full set of a^ m 's and the spectral coefficients Cg. We show that using 
an Equidistant Cylindrical Projection of the sky substantially speeds up the calculations. 
Thus, generating and/or inverting a full sky, high resolution map can be easily achieved 
with present day computer technology. 



Subject Headings: cosmic microwave background 



1 Introduction 



The angular power spectrum of CMB anisotropies is a gold-mine of cosmological informa- 
tion. It sensitively depends upon a number of parameters: the total and baryonic density 
parameters, Q an d &>b', the cosmological constant, A; the Hubble constant, H ; the spectral 
indices, n s and n t , and the amplitudes of scalar and tensor metric fluctuations; the redshift, 
z rh , at which the universe could have been reionized. Because of their planned high sensi- 
tivity and high angular resolution, future space missions will measure the anisotropy power 
spectrum with great accuracy. Thus, all these cosmological parameters will be determined 
with an unprecedented precision (Bersanelli etal. 1996; Jungman etal. 1996) 

To achieve these goals, Monte Carlo simulations of the CMB anisotropy pattern have 
been and will be more and more crucial in this game. From one hand, they allow to 
prepare a mission, to optimize the observational strategy and to test for different payload 
configurations. On the other hand, they are important for the data analysis and, for 
example, to look for systematics. 

Up to now, these simulations were realized without problems. Experiments with high 
angular resolution observed only very limited region of the sky: FFT techniques easily 
provided several realizations of small (and, hence, flat) patches of the sky (see e.g. Kogut, 
Hinshaw and Bennett, 1995). Experiments with large sky coverage, such as COBE/DMR, 
had low resolution: several full sky, CMB anisotropy maps were easily generated through 
a spherical harmonic expansion with a low ( ^ 100) number of harmonics. 

A potential problem it is claimed to arise in generating even a single high resolution, 
full sky map (see, e.g. Saez, Holtmann and Smoot, 1996), too much, it is believed, for 
present computer technology. It is generally perceived as a heavy, almost impossible, 
computational task also the inverse problem, that is extracting out of an observed full sky, 
high resolution map the anisotropy spectrum up to £ ^ 1000. In fact, to our knowledge, 
anisotropy spectrum estimates have always been done applying FFT techniques either to 
small patches of , or even to the whole (FFT simulated) sky. The latter approach is in 
principle wrong and its accuracy has still to be checked for. 
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The purpose of this paper is quite technical, but of interest, we believe, to the large 
community involved in future CMB anisotropy experiments. We want to present a fast 
algorithm for: i) generating high resolution ( ^ 10'), full sky maps; ii) reconstructing all 
the coefficients of a spherical harmonic expansion, and hence the spectral coefficient Ce, up 
to £ ^ 1000. Both these tasks can be easily achieved on currently available workstations. 
Thus, the plane of this Letter is as follows. In Sect. 2 we will describe the method. In 
Sect. 3 we will discuss numerical results and the efficiency of the method. Finally in Sect. 4 
we will present a brief summary of our main findings. 

2 Method 

Generating a CMB anisotropy map is in principle very simple. The temperature fluctuation 
observed along a line of sight, 7, can be conveniently described by a spherical harmonic 
expansions: 

imax i 

— (X^)=Y, J2 a £m(x) Y £m(l) (!) 
1 £=0 ra=-l 

where a^_ m = (— l) m a* lm . The ai m 's are random variables of the observer position, x, 
gaussian distributed (at least in most of the inflation based scenarios), with zero mean and 
variances (\ae m \ 2 ) = Cg. In simulating the CMB, primary anisotropy pattern the sum over 
£ usually starts from two. In fact, from one hand the monopole vanishes by construction, 
being the mean (over the sky) CMB anisotropy. On the other hand, the dipole components 
are dominated by the Doppler anisotropy, induced by our peculiar motion relative to the 
comoving frame (see e.g. Kogut etal. 1994). On the same line, we will keep the sum over 
£ from to £ ma x, but we set a 00 = a\^\ = a^o = = 0. 

The CV's are the main prediction of a theory of structure formation (see e.g. Hu and 
Sugiyama, 1996). Thus, for a given scenario {i.e. for given CVs) and for a given statistics, 
we can generate a random set of a£ m 's and, from Eq.([l|), a CMB anisotropy map. In 
practice, using the spherical harmonic expansion as in Eq.(|I|) is not very efficient: for each 
line of sight 7 we should evaluate Ye m (^) for each value of £ and m. Fortunately, it is 
possible to rewrite Eq.([TJ) in a form more suitable for numerical implementation. 
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First, it is easy to verify that the double sum in Eq.(|l|), 



l=lmax -<r^m=+£ 
1=Iq 2-im=—ti 



is completely 



equivalent to 



m=+£r, 
m=—t v 



ax Y^p=\Z\ x '- we sample exactly the same region of the £ — m space, by 



columns (in the former) or by rows (in the latter case; c.f. Fig.l). Second, let us write 



f.m\ 



\7(cos6)e 



im<t> 



(2) 



where 



a; 



2£+l (£-m)\ 



PI" (cos f 



(3) 



\\ 4tt (£ + m)\ 

\J m = (— ) m A™ and P™(cos9) are the associated Legendre polynomials. It is then possible 
to rewrite Eq.(P as follows: 



AT 

-r 1 



^ b m (0)e im4 



where 



t=\m\ 



\ m 



(4) 



(5) 



and 6_ m = b* m . Written in this way, Eq.(|) highlights a couple of attractive features. If 
we use an Equidistant Cylindrical Projection (hereafter ECP) of the sky, which conserves 
distances along meridians and along the equator, the anisotropy map can be thought 
as a rectangular matrix of N$ times Ng(= N^/2) squared pixels, each of dimension ~ 
20' x 20'(1024/A' r ^) 2 . In this projection, the temperature anisotropy along parallels (i.e. 
at fixed 9) is nothing more than the ID Fourier transform of the coefficients b m (9) 1 s [c.f. 
Eq.(|])], very efficiently computed with FFT techniques. If we regard the sum of Eq.(^j) as 
a Fourier expansion, then t max must be fixed to be as it plays the role of the Nyquist 

critical frequency of the problem. Second, the A™ are evaluated by standard recurrence 
relations: 



a: 



\ m 
A m+1 



1 2m +1 (2m- V 



A-n 

xV2m + 3A™ 



(2m)! 



[l-x 



2\m/2 



(l + m- I) (I - m - 1) 
V (2£-3)(2£- 1) 



\ m 
A t-2 



U 2 



£ 2 — m 2 



(6) 



where x = cos 9. Because of these relations, the b m (9) 's can be computed very efficiently, 
as we can perform the sum in Eq.([5]) while computing the A™'s. Such a computation is 
further simplified because \i m (cos9) = ±A^ m [cos(7r — 9)], the plus (minus) sign holding if 
I and m are (are not) both even or both odd. 

Generating the 6 m (#)'s at fixed 9 requires evaluating « (oc iV|) recurrence relations 
{imax — m recurrence relations for each value of m). The CPU time needed for FFT-ing 
the b m (9) , s in principle scales as N^lnN^. However, the FFT is so fast that for < 4096 
most of the time is spent for generating the b m {9) 's. Finally, we have to evaluate the 
b m (9)'s Ng times. So, at the end, the total CPU time needed for generating an ECP of the 
anisotropy pattern is expected to scale as iV|. 

At this point it is quite easy to address also the inverse problem, that of evaluating from 
an observed high resolution, full sky map the set of coefficients ae m . It is very well known 
that the orthonormality of the spherical harmonics allows to invert Eq.([I]) and write: 

r AT 

ai m = J dn—(j)Y; m (j) (7) 

This sounds awful: in principle for each i and m we should evaluate Y^ m (7) for a given 7 
and integrate over the whole sky. Fortunately, after substituting Eq.(^) in Eq.(|7]) we can 
write: 

a lm = J sin 9d9 XT (9)b rn (9) (8) 

where 

b m (9) = d(pA(4>,9)exp(-imcf>) (9) 
Jo 

Thus, Eq.(P) is the conjugate of Eq.(|J): the 6 m 's are the Fourier anti-transform of the 
anisotropy pattern along a parallel in the ECP of the sky, and are easily computed at fixed 
9 with a FFT. In conclusion, inverting a map to obtain the a£ m 's requires ~ £ max (oc NT) 
recurrence relations for evaluating the A™'s, plus a FFT to evaluate the 6 m 's. All this 
must be done Ng (oc N^) times to be able to perform the integral in Eq.(^). Using these 
tricks, we can invert a full sky, high resolution map with CPU times which scale as iV| 
(the evaluation of the 6^ m 's is basically instantaneous), and in principle comparable with 
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those needed for generating a map. As in that case, the actual CPU time can be further 
reduced by exploiting the simmetries of the A^ m 's evaluated at 6 and -n — 9, respectively. 

3 Numerical Results 

In the previous Section we described an algorithm for generating and/or inverting a high 
resolution, full sky map of the CMB anisotropy. In this Section we will briefly discuss the 
actual performances of this algorithm on a DEC 1000/200. We will consider hereafter the 
standard Cold Dark Matter model. 

In Fig. 2 we plot the CPU time needed for generating a full sky, ECP of the CMB 
anisotropy pattern as a function of the angular resolution, the only free parameter we can 
play with. In fact, we sample the ECP of the sky with x N d squared pixels of dimension 
20' x 2O / (lO24/jV ) 2 . This fixes N e = t max = The expected scaling with JVj (see 

Sect. 2) is recovered with good precision. We want to stress that only lh of CPU time is 
needed to generate a full sky, ECP of the CMB anisotropy with a resolution of 5' (i.e. 
N<f) = 4096). 

Using an ECP is not a limitation. In fact, once we obtain an ECP of the anisotropy 
pattern, we can reproduce it in any given projection. As an example we show in Plate 1 an 
ECP, 10' resolution map (obtained in only 8 minutes of CPU time) and the corresponding 
Equal Area Projection (hereafter EAP). The latter is obtained by the former using standard 
spherical trigonometry. We verified that this procedure is numerically stable. In fact, if we 
transform an ECP to an EAP and the obtained EAP back to an ECP, we reproduce the 
initial anisotropy pattern exactly. Thus, from observations of the CMB anisotropy we can 
create an ECP of the sky and then apply our inversion algorithm. In Fig. 2 we also show 
the CPU time needed for inverting an ECP map as a function of the map resolution. The 
CPU time scales roughly as iV|. Again, ~ lh of CPU time is needed to recover from a 5' 
resolution map the entire set of a^ m 's, roughly 4 • 10 6 coefficients. 

In Figs. 3 and 4 we show the precision of our inversion algorithm for £ max = 1024 
(corresponding to a resolution of 10') . The percentage error between the recovered a^ m 's 
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and the input ones is large only for £ ~ £ ma x and m ~ few, a very small portion of the 
allowed region of the i — m space. This is due to the fact that for m — > 0, A° oc Pi(cos6): 
for large values of £, this is a highly oscillating function of 9. On the contrary, for m — > £, 
\\ oc (1 — cos#)^ 2 , a quite smooth function of the azimuthal angle [c.f. Eq. ([])]. So, a 
simple trapezoidal rule for performing the integral along meridians in Eq.® gives a poor 
result only for very large values of I and quite small values of m. However, this is not a 
crucial problem. In fact, we are mostly interested in evaluating the spectral coefficients Cf . 
These are obtained from the recovered ci£ m 's as follows: 

1 1 

/^estimated \ * i |2 /ia\ 

1 = 2lTi 1 - P laeml (10) 

m=—t 

It is clear that the error we make in recovering the a^ m 's for large i and small m is 
highly diluted in the sum of Eq. fllOl) . In Fig. 5 we show the (percentage) error between 
the recovered and the input C/s as a function, again, of the map resolution. The recovered 
spectrum has a maximum error of ~ 0.1% up to I ^ 1500 for = 4096, i.e. for a pixel 
size of 5' x 5'. 

4 Conclusions 

We present a fast algorithm for: i) generating high resolution, full sky maps of the CMB 
anisotropy; ii) evaluating out of an observed map the a^ m 's and then the spectral coefficient 
Cg's. The basic trick for speeding up the calculation consists in generating and/or inverting 
a full sky map using an ECP. It is this projection that allows the use of a FFT either in 
Eq.(f|) and/or in Eq.(|9]). If we are interested in probing the anisotropy power spectrum 
up to £ max ~ 1000, then = 2£ max = 2048 and we need only 8 minutes of CPU time 
for either generating or inverting a 10' resolution, full sky map. Pushing the sampling 
down to 5' boosts the needed CPU time up to one hour. Our algorithm also allows us to 
fully exploit a parallel architecture, such as the one of APEmille (Bartoloni etal. 1995). 
This machine is composed by 1024 processors, each of them slower by roughly a factor 
of two w.r.t. a DEC 1000/200. So, with such a machine one should be able to produce 
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1024, 5' resolution maps of the anisotropy pattern in a couple of hours. Details about this 
application will be discussed elsewhere. 

In addressing the problem of inverting a full sky map, we assumed full sky coverage and 
we fully exploited the orthonormality of the spherical harmonics. We test our algorithm 
against a pure CMB anisotropy pattern. In the realistic case, a /z-wave map will be 
the superposition of different processes (CMB anisotropy, Galactic foregrounds, secondary 
anisotropy due to clusters of galaxies, point sources, etc.) and the sky coverage can be not 
complete. The separation of the CMB anisotropy pattern from Galactic and extragalactic 
foregrounds has been studied in details (Bouchet etal. 1994; Bouchet etal. 1995; Bersanelli 
etal. 1996), but considering only small (10° x 10°) patches of the sky. We will discuss an 
application of our algorithm to the problem of foreground subtraction and not complete 
sky coverage in a forthcoming paper. 
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Figure 1: The filled squares show (for £ max — 3) the portion of the £ — m space probed 
by the a^ m 's. We can recover all the allowed £ — m pairs either by moving along columns 

l^=o i^ m =-e) 01 b Y rows {l^ m =-e max 2^e=\m\ )• 
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Figure 2: The CPU time needed, on a DEC 1000, to generate a full sky map (solid line) or 
to invert a map to generate the full set of a^ m 's (dotted line) as a function of the sampling 
of an Equidistant Cylindrical Projection of the sky. 
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Figure 3: The percentage error between the recovered and the input, real part of the ag m 's. 
The labels indicate the 0.01% and the 0.1% isolevel, respectively. Only in a very small 
portion of the I — m plane the error is larger than 10% (dashed line). 
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Figure 4: Same as in Fig.3, but for the imaginary part of the a£ m 's. 
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Figure 5: The percentage error between the recovered and the input spectral coefficients 
C/s as a function of £ for different resolutions of the Equidistant Cylindrical Projection 
of the CMB anisotropy pattern. Dotted, dashed, dot-dashed and continous lines refer to 
map resolution of ~ 40' (£ max = 256), ~ 20' (£ max = 512), ~ 10' (£ max = 1024), ~ 5' 
{£max = 2048), respectively. 



Figure 6: A realization of the CMB anisotropy pattern in a standard Cold Dark Matter 
model in the Equidistant Cylindrical (panel a) and in the Equal Area (panel b) Projections. 
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This figure "fig6.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/astro-ph/9703084vl 



